
## Install & load packages (all at once)
list.of.packages <- c("classInt","flexmix","MCMCpack","BMA","mvtnorm","maptools","spdep","rgeos","mgcv","mvtnorm","sfsmisc","xtable")
new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]; if(length(new.packages)){install.packages(new.packages,dependencies=TRUE)}
lapply(list.of.packages, require, character.only = TRUE)


################################
# Interaction BMA
################################

bma.int <- function(mod,ix){

int. <- c(); nint. <- c()
for(k in 1:length(ix)){
int. <- c(int.,which(grepl(ix[[k]][1],mod$label)&grepl(ix[[k]][2],mod$label)&grepl(ix[[k]][2],mod$label)))
nint. <- c(nint.,which(grepl(ix[[k]][1],mod$label)&(!(grepl(ix[[k]][2],mod$label)&grepl(ix[[k]][3],mod$label)))))
}
int. <- sort(unique(int.))
nint. <- sort(unique(nint.))

# Old QOI's
postprob.old <- mod$postprob
probne0.old <- apply(mod$which,2,function(x){weighted.mean(x,mod$postprob)})
postmean.old <- apply(mod$mle,2,function(x){weighted.mean(x,mod$postprob)})
postse.old <- sqrt(abs(apply(mod$se^2+mod$mle^2,2,function(x){weighted.mean(x,mod$postprob)})-postmean.old^2))

# New QOI's
postprob.new <- mod$postprob[-nint.]/sum(mod$postprob[-nint.])
sum(postprob.new)
probne0.new <- apply(mod$which[-nint.,],2,function(x){weighted.mean(x,postprob.new)})
postmean.new <- apply(mod$mle[-nint.,],2,function(x){weighted.mean(x,postprob.new)})
postse.new <- sqrt(abs(apply(mod$se[-nint.,]^2+mod$mle[-nint.,]^2,2,function(x){weighted.mean(x,postprob.new)})-postmean.new^2))

new.mod <- list(namesx=mod$namesx,probne0=probne0.new*100,postmean=postmean.new,postsd=postse.new)

return(new.mod)
}


#ix <- list(c("MCH_RUS","IND_MCH","RUS_MAJ"),c("MIN_RUS","IND_MIN","RUS_MAJ"),c("MTL_RUS","IND_MTL","RUS_MAJ"))
#new.mod <- bma.int(mod=mod,ix=ix)


#############################
## BMA plot                ##
#############################




BMA.plot <- function(mod,labz=names(mod$probne0),dvlab=NULL,dv.cex=1.5,n=5000){
set.seed(1)
betas <- rmvnorm(n,mod$postmean[-1],diag(mod$postsd[-1]^2))
par(mar=c(6,9,ifelse(length(dvlab)>0,3,.5),0.5),xpd=F)
plot(x=0,y=0,col=NA,xlim=c(floor(min(betas)),ceiling(max(betas))),ylim=c(0,ncol(betas)+1),xaxt="n",yaxt="n",ylab="",xlab="",bty="n")
eaxis(1, padj=-0.5, cex.axis=0.8)
rect(xleft=floor(min(betas)), ybottom=0, xright=ceiling(max(betas)), ytop=ncol(betas)+1,col="grey90",border=NA)
segments(x0=floor(min(betas)):ceiling(max(betas)),x1=floor(min(betas)):ceiling(max(betas)),y0=-1,y1=ncol(betas)+1,lty=3,col="darkgrey")
for(i in 1:ncol(betas)){points(x=betas[,ncol(betas)-i+1],y=runif(n,i-.4,i+.4),col=rgb(.1,.3,.5,mod$probne0[ncol(betas)-i+1]*.005),pch=16,cex=.05)
segments(x0=mean(betas[,ncol(betas)-i+1]),x1=mean(betas[,ncol(betas)-i+1]),y0=i-.5,y1=i+.5,col=rgb(1,1,1),lwd=2)
segments(x0=quantile(betas[,ncol(betas)-i+1],.025),x1=quantile(betas[,ncol(betas)-i+1],.025),y0=i-.5,y1=i+.5,lty=1,col=rgb(1,1,1))
segments(x0=quantile(betas[,ncol(betas)-i+1],.975),x1=quantile(betas[,ncol(betas)-i+1],.975),y0=i-.5,y1=i+.5,lty=1,col=rgb(1,1,1))
}
segments(x0=0,x1=0,y0=-1,y1=ncol(betas)+1,col="red",lwd=2,lty=1)
axis(2,ncol(betas):1,labels=labz,las=1,cex.axis=.8)
par(xpd=NA)
cutz <- seq(floor(min(betas)),ceiling(max(betas)),length.out=11)
j <- 1
for(j in 1:10){
rect(xleft=cutz[j], ybottom=-ncol(betas)/4, xright=cutz[j+1], ytop=-ncol(betas)/5,col=rgb(.1,.3,.5,c(1:10/10)[j]),border=NA,cex=.8)
}
text(x=c(floor(min(betas)),mean(c(floor(min(betas)),ceiling(max(betas)))),ceiling(max(betas))),y=-ncol(betas)/3.5,labels=c(0,.5,1),cex=.8)
axis(2,at=-ncol(betas)/4,labels=expression(P(beta!=0)),las=1,tick=F,cex.axis=.8)
if(length(dvlab)>0){
axis(3,at=mean(c(floor(min(betas)),ceiling(max(betas)))),labels=dvlab,tick=F,cex.axis=dv.cex,padj=.5)
}
par(xpd=F)
}

#########################

BMA.plot.surv <- function(mod,labz=mod$namesx,dvlab=NULL,dv.cex=1.5,n=5000,xlim=NULL){
set.seed(1)
betas <- exp(rmvnorm(n,mod$postmean,diag(mod$postsd^2)))
if(length(xlim)==0){xlim<-c(floor(min(betas))+.001,ceiling(max(betas)))}
#if(length(xlim)==0){xlim<-c(floor(min(betas)*10)/10+.001,ceiling(max(betas)))}
par(mar=c(6,9,ifelse(length(dvlab)>0,3,.5),0.5),xpd=F)
plot(x=0,y=0,col=NA,xlim=xlim,ylim=c(0,ncol(betas)+1),yaxt="n",xaxt="n",ylab="",xlab="",log="x",cex.axis=.8,bty="n")
eaxis(1, padj=-0.5, cex.axis=0.8)
rect(xleft=floor(min(betas))+.001, ybottom=0, xright=ceiling(max(betas)), ytop=ncol(betas)+1,col="grey90",border=NA)
segments(x0=10^(floor(log(floor(min(betas))+.001,10)):ceiling(log(floor(max(betas))+.001,10))),x1=10^(floor(log(floor(min(betas))+.001,10)):ceiling(log(floor(max(betas))+.001,10))),y0=-1,y1=ncol(betas)+1,lty=3,col="darkgrey")
for(i in 1:ncol(betas)){points(x=betas[,ncol(betas)-i+1],y=runif(n,i-.4,i+.4),col=rgb(.1,.3,.5,mod$probne0[ncol(betas)-i+1]*.005),pch=16,cex=.05)
segments(x0=mean(betas[,ncol(betas)-i+1]),x1=mean(betas[,ncol(betas)-i+1]),y0=i-.5,y1=i+.5,col=rgb(1,1,1),lwd=2)
segments(x0=quantile(betas[,ncol(betas)-i+1],.025),x1=quantile(betas[,ncol(betas)-i+1],.025),y0=i-.5,y1=i+.5,lty=1,col=rgb(1,1,1))
segments(x0=quantile(betas[,ncol(betas)-i+1],.975),x1=quantile(betas[,ncol(betas)-i+1],.975),y0=i-.5,y1=i+.5,lty=1,col=rgb(1,1,1))
}
segments(x0=1,x1=1,y0=-1,y1=ncol(betas)+1,col="red",lwd=2,lty=1)
axis(2,ncol(betas):1,labels=labz,las=1,cex.axis=.8)
par(xpd=NA)
cutz <- exp(seq(log(floor(min(betas))+.001),log(ceiling(max(betas))),length.out=11))
j <- 1
for(j in 1:10){
rect(xleft=cutz[j], ybottom=-ncol(betas)/4, xright=cutz[j+1], ytop=-ncol(betas)/5,col=rgb(.1,.3,.5,c(1:10/10)[j]),border=NA)
}
text(x=exp(c(log(floor(min(betas))+.001),mean(c(log(floor(min(betas))+.001),log(ceiling(max(betas))))),log(ceiling(max(betas))))),y=-ncol(betas)/3.5,labels=c(0,.5,1),cex=.8)
axis(2,at=-ncol(betas)/4,labels=expression(P(exp(beta)!=1)),las=1,tick=F,cex.axis=.8)
if(length(dvlab)>0){
axis(3,at=exp(mean(c(log(floor(min(betas))+.001),log(ceiling(max(betas)))))),labels=dvlab,tick=F,cex.axis=dv.cex,padj=.5)
}
par(xpd=F)
}




BMA.plot.surv2 <- function(mod,labz=mod$namesx,dvlab=NULL,dv.cex=1.5,n=5000,xlim=NULL){
set.seed(1)
betas <- exp(rmvnorm(n,mod$postmean,diag(mod$postsd^2)))
#if(length(xlim)==0){xlim<-c(floor(min(betas))+.001,ceiling(max(betas)))}
if(length(xlim)==0){xlim<-c(floor(min(betas)*10)/10+.001,ceiling(max(betas)))}
par(mar=c(6,9,ifelse(length(dvlab)>0,3,.5),0.5),xpd=F)
plot(x=0,y=0,col=NA,xlim=xlim,ylim=c(0,ncol(betas)+1),yaxt="n",xaxt="n",ylab="",xlab="",log="x",cex.axis=.8,bty="n")
eaxis(1, padj=-0.5, cex.axis=0.8)
rect(xleft=floor(min(betas))+.001, ybottom=0, xright=ceiling(max(betas)), ytop=ncol(betas)+1,col="grey90",border=NA)
segments(x0=10^(floor(log(floor(min(betas))+.001,10)):ceiling(log(floor(max(betas))+.001,10))),x1=10^(floor(log(floor(min(betas))+.001,10)):ceiling(log(floor(max(betas))+.001,10))),y0=-1,y1=ncol(betas)+1,lty=3,col="darkgrey")
for(i in 1:ncol(betas)){points(x=betas[,ncol(betas)-i+1],y=runif(n,i-.4,i+.4),col=rgb(.1,.3,.5,mod$probne0[ncol(betas)-i+1]*.005),pch=16,cex=.05)
segments(x0=mean(betas[,ncol(betas)-i+1]),x1=mean(betas[,ncol(betas)-i+1]),y0=i-.5,y1=i+.5,col=rgb(1,1,1),lwd=2)
segments(x0=quantile(betas[,ncol(betas)-i+1],.025),x1=quantile(betas[,ncol(betas)-i+1],.025),y0=i-.5,y1=i+.5,lty=1,col=rgb(1,1,1))
segments(x0=quantile(betas[,ncol(betas)-i+1],.975),x1=quantile(betas[,ncol(betas)-i+1],.975),y0=i-.5,y1=i+.5,lty=1,col=rgb(1,1,1))
}
segments(x0=1,x1=1,y0=-1,y1=ncol(betas)+1,col="red",lwd=2,lty=1)
axis(2,ncol(betas):1,labels=labz,las=1,cex.axis=.8)
par(xpd=NA)
#cutz <- exp(seq(log(floor(min(betas))+.001),log(ceiling(max(betas))),length.out=11))
cutz <- exp(seq(log(xlim[1]),log(xlim[2]),length.out=11))
j <- 1
for(j in 1:10){
rect(xleft=cutz[j], ybottom=-ncol(betas)/4, xright=cutz[j+1], ytop=-ncol(betas)/5,col=rgb(.1,.3,.5,c(1:10/10)[j]),border=NA)
}
text(x=exp(c(log(floor(min(betas)*10)/10+.001),mean(c(log(floor(min(betas)*10)/10+.001),log(ceiling(max(betas))))),log(ceiling(max(betas))))),y=-ncol(betas)/3.5,labels=c(0,.5,1),cex=.8)
axis(2,at=-ncol(betas)/4,labels=expression(P(exp(beta)!=1)),las=1,tick=F,cex.axis=.8)
if(length(dvlab)>0){
axis(3,at=exp(mean(c(log(floor(min(betas)*10)/10+.001),log(ceiling(max(betas)))))),labels=dvlab,tick=F,cex.axis=dv.cex,padj=.5)
}
par(xpd=F)
}



#########################
BMA.mplot <- function(mod,n.sim=1000,n.out=100,rang=c(min(mod$x[,var]),max(mod$x[,var])),var=NULL,xlab=NULL,ylab=NULL,rugz=FALSE,barz=TRUE,sq=NULL){
X.sim <- as.data.frame(matrix(apply(mod$x,2,median),nrow=n.out,ncol=ncol(X),byrow=T))
names(X.sim) <- names(mod$x)
X.sim[,var] <- seq(rang[1],rang[2],length.out=n.out)
if(length(sq)>0){X.sim[,sq] <- X.sim[,var]^2}
set.seed(1)
betas <- rmvnorm(n.sim,mod$postmean[-1],diag(mod$postsd[-1]^2))
pmat <- matrix(NA,nrow=n.out,ncol=n.sim)
for(i in 1:n.sim){
pmat[,i] <- mod$linkinv(as.matrix(X.sim)%*%betas[i,])
}
y.mean <- apply(pmat,1,function(x){mean(x)})
y.lo <- apply(pmat,1,function(x){quantile(x,.025)})
y.hi <- apply(pmat,1,function(x){quantile(x,.975)})

if(length(xlab)==0){xlab=var}
if(length(ylab)==0){ylab="Outcome"}

par(mar=c(4,4,0.5,0.5))
plot(x=0,y=0,col=NA,ylim=c(min(y.lo),max(y.hi)),xlim=c(min(X.sim[,var]),max(X.sim[,var])),bty="l",xlab=xlab,ylab=ylab,cex.lab=.8)
rect(xleft=min(X.sim[,var]), ybottom=min(y.lo), xright=max(X.sim[,var]), ytop=max(y.hi),col="grey90",border=NA)
if(barz==TRUE){
segments(x0=X.sim[,var],x1=X.sim[,var],y0=y.lo,y1=y.hi,col=rgb(.1,.3,.5,.75))
}
if(barz==FALSE){
polygon(x=X.sim[c(1:n.out,n.out:1),var],y=c(y.lo,rev(y.hi)),col=rgb(.1,.3,.5,.75),border=NA)
}
lines(x=X.sim[,var],y=y.mean,lwd=2)
if(rugz==TRUE){rug(jitter(mod$x[,var]))}
}




#########################
BMA.mplot2 <- function(mod,n.sim=1000,n.out=100,rang=c(min(mod$x[,var]),max(mod$x[,var])),var=NULL,xlab=NULL,ylab=NULL,rugz=FALSE,barz=TRUE,sq=NULL,cexz=.8,xpos="topleft",ypos=NULL){
  X.sim <- as.data.frame(matrix(apply(mod$x,2,median),nrow=n.out,ncol=ncol(X),byrow=T))
  names(X.sim) <- names(mod$x)
  X.sim[,var] <- seq(rang[1],rang[2],length.out=n.out)
  if(length(sq)>0){X.sim[,sq] <- X.sim[,var]^2}
  set.seed(1)
  betas <- rmvnorm(n.sim,mod$postmean[-1],diag(mod$postsd[-1]^2))
  pmat <- matrix(NA,nrow=n.out,ncol=n.sim)
  for(i in 1:n.sim){
    pmat[,i] <- mod$linkinv(as.matrix(X.sim)%*%betas[i,])
  }
  y.mean <- apply(pmat,1,function(x){mean(x)})
  y.lo <- apply(pmat,1,function(x){quantile(x,.025)})
  y.hi <- apply(pmat,1,function(x){quantile(x,.975)})
  
  if(length(xlab)==0){xlab=var}
  if(length(ylab)==0){ylab="Outcome"}
  
  par(mar=c(3,2,0.5,0.5))
  plot(x=0,y=0,col=NA,ylim=c(min(y.lo),max(y.hi)),xlim=c(min(X.sim[,var]),max(X.sim[,var])),bty="l",xlab="",ylab=ylab,cex.lab=cexz)
  title(xlab = xlab, cex.lab = cexz,line = 2)
  rect(xleft=min(X.sim[,var]), ybottom=min(y.lo), xright=max(X.sim[,var]), ytop=max(y.hi),col="grey90",border=NA)
  if(barz==TRUE){
    segments(x0=X.sim[,var],x1=X.sim[,var],y0=y.lo,y1=y.hi,col=rgb(.1,.3,.5,.75),lwd=1.5)
  }
  if(barz==FALSE){
    polygon(x=X.sim[c(1:n.out,n.out:1),var],y=c(y.lo,rev(y.hi)),col=rgb(.1,.3,.5,.75),border=NA)
  }
  lines(x=X.sim[,var],y=y.mean,lwd=2)
  if(rugz==TRUE){rug(jitter(mod$x[,var]))}
  legend(x=xpos,y=ypos,title=ylab,fill=rgb(.1,.3,.5,.75),legend="95% credible interval",bty="n",cex=cexz)
}



#########################
BMA.pred <- function(mod,n.sim=1000,n.out=100,rang=c(min(mod$x[,var]),max(mod$x[,var])),var=NULL,sq=NULL){
X.sim <- as.data.frame(matrix(apply(mod$x,2,median),nrow=n.out,ncol=ncol(X),byrow=T))
names(X.sim) <- names(mod$x)
X.sim[,var] <- seq(rang[1],rang[2],length.out=n.out)
if(length(sq)>0){X.sim[,sq] <- X.sim[,var]^2}
set.seed(1)
betas <- rmvnorm(n.sim,mod$postmean[-1],diag(mod$postsd[-1]^2))
pmat <- matrix(NA,nrow=n.out,ncol=n.sim)
for(i in 1:n.sim){
pmat[,i] <- mod$linkinv(as.matrix(X.sim)%*%betas[i,])
}
y.mean <- apply(pmat,1,function(x){mean(x)})
y.lo <- apply(pmat,1,function(x){quantile(x,.025)})
y.hi <- apply(pmat,1,function(x){quantile(x,.975)})
return(list(y.mean=y.mean,y.lo=y.lo,y.hi=y.hi))
}



#########################


BMA.X <- function(X,xvals,varnames,ix=NULL){
X.sim <- as.data.frame(matrix(apply(X,2,median),nrow=1,ncol=ncol(X),byrow=T))
names(X.sim) <- varnames
for(k in 1:length(xvals)){
X.sim[,xvals[[k]][1]] <- as.numeric(xvals[[k]][2])
}
if(length(ix)>0){
for(k in 1:length(ix)){
X.sim[,ix[[k]][1]] <- X.sim[,ix[[k]][2]]*X.sim[,ix[[k]][3]]
}}
return(X.sim)
}


#######################


BMA.pred.HR <- function(mod,n.sim=1000,X.0=NULL,X.1=NULL){
X.sim <- X.1-X.0
#set.seed(1)
betas <- rmvnorm(n.sim,mod$postmean,diag(mod$postsd^2))
pmat <- matrix(NA,nrow=1,ncol=n.sim)
for(i in 1:n.sim){
pmat[,i] <- as.matrix(X.sim)%*%betas[i,]
}
hr.mean <- exp(median(pmat))
hr.lo <- exp(quantile(pmat,.025))
hr.hi <- exp(quantile(pmat,.975))
#pmat <- exp(pmat)
#hr.mean <- median(pmat)
#hr.lo <- quantile(pmat,.025)
#hr.hi <- quantile(pmat,.975)
return(list(hr.mean=hr.mean,hr.lo=hr.lo,hr.hi=hr.hi))
}



#######################


BMA.pred.HR2 <- function(mod,n.sim=1000,X.0=NULL,X.1=NULL,rnd=2){
X.sim <- X.1-X.0
#set.seed(12345)
betas <- rmvnorm(n.sim,mod$postmean,diag(mod$postsd^2))
pmat <- matrix(NA,nrow=1,ncol=n.sim)
for(i in 1:n.sim){
pmat[,i] <- as.matrix(X.sim)%*%betas[i,]
}
hr.mean <- (exp(median(pmat))-1)*100
hr.lo <- (exp(quantile(pmat,.025))-1)*100
hr.hi <- (exp(quantile(pmat,.975))-1)*100
#hr.mean <- (median(exp(pmat))-1)*100
#hr.lo <- (quantile(exp(pmat),.025)-1)*100
#hr.hi <- (quantile(exp(pmat),.975)-1)*100
return(list(hr.mean=round(hr.mean,rnd),hr.ci=paste0("(",round(hr.lo,rnd),", ",round(hr.hi,rnd),")")))
}



#######################



BMA.pred.alt <- function(mod,n.sim=1000,X.0=NULL,X.1=NULL,invlink,rnd=2){
#set.seed(1)
betas <- rmvnorm(n.sim,mod$postmean[-1],diag(mod$postsd[-1]^2))
pmat1 <- matrix(NA,nrow=1,ncol=n.sim)
pmat0 <- matrix(NA,nrow=1,ncol=n.sim)
for(i in 1:n.sim){
pmat1[,i] <- invlink(as.matrix(X.1)%*%betas[i,])
pmat0[,i] <- invlink(as.matrix(X.0)%*%betas[i,])
}
y0 <- c(round(median(pmat0),rnd),paste0("(",round(quantile(pmat0,.025),rnd),",",round(quantile(pmat0,.975),rnd),")"))
y1 <- c(round(median(pmat1),rnd),paste0("(",round(quantile(pmat1,.025),rnd),",",round(quantile(pmat1,.975),rnd),")"))
ate <- c(round(median(pmat1-pmat0),rnd),paste0("(",round(quantile(pmat1-pmat0,.025),rnd),",",round(quantile(pmat1-pmat0,.975),rnd),")"))
rr <- c(round(median((pmat1/pmat0-1)*100),rnd),paste0("(",round(quantile((pmat1/pmat0-1)*100,.025),rnd),",",round(quantile((pmat1/pmat0-1)*100,.975),rnd),")"))
out <- as.data.frame(rbind(y0,y1,ate,rr))
names(out) <- c("Mean","CI")
return(out)
}



#############################
## Custom palette function ## 
#############################

custom.palette <- function(var,col.vec,n.color=5,choose=F,style="pretty",fixedBreaks){
	
col.ramp <- colorRampPalette(col.vec, space = "rgb")


# Create palette with n color classes
pal <- col.ramp(n=n.color) 

## (1) Show options
if(choose==T){
#classes_fx <- classIntervals(var, n=n.color, style="fixed", fixedBreaks=pretty(var,min.n=3,n=n.color), rtimes = 1)
classes_pr<-classIntervals(var,n=n.color,style="pretty",rtimes=1)
classes_sd<-classIntervals(var,n=n.color,style="sd",rtimes=1)
classes_fi<-classIntervals(var,n=n.color,style="fisher",rtimes=3)
classes_eq<-classIntervals(var,n=n.color,style="equal",rtimes=1)
classes_km<-classIntervals(var,n=n.color,style="kmeans",rtimes=1)
classes_qt<-classIntervals(var,n=n.color,style="quantile",rtimes=1)
classes_hc<-classIntervals(var,n=n.color,style="hclust",rtimes=1)
classes_bc<-classIntervals(var,n=n.color,style="bclust",rtimes=1)

par(mar=c(2,2,2,1)+0.1, mfrow=c(2,4))
#plot(classes_fx, pal=pal, main="Fixed Intervals", xlab="", ylab="")
plot(classes_pr, pal=pal, main="pretty", xlab="", ylab="")
plot(classes_sd, pal=pal, main="sd", xlab="", ylab="")
plot(classes_fi, pal=pal, main="fisher (jenks)", xlab="", ylab="")
plot(classes_km, pal=pal, main="kmeans", xlab="", ylab="")
plot(classes_eq, pal=pal, main="equal", xlab="", ylab="")
plot(classes_qt, pal=pal, main="quantile", xlab="", ylab="")
plot(classes_hc, pal=pal, main="hclust", xlab="", ylab="")
plot(classes_bc, pal=pal, main="bclust", xlab="", ylab="")
}

## (2) Create cupoints
if(choose==F){
	
## Define intervals
if(style!="fixed"){classes<-classIntervals(var,n=n.color,style=style,rtimes=1)}
if(style=="fixed"){classes <- classIntervals(var, n=n.color, style="fixed", fixedBreaks=fixedBreaks, rtimes = 1)}

## Create vector of colors
cols <- findColours(classes, pal)
return(cols)	
}
}


####################################
## round.intervals                ## 
####################################


round.intervals <- function(cols,rnd){
labs <- names(attr(cols, "table"))
labs <- gsub("[","",labs,perl=F,fixed=T)
labs <- gsub("]","",labs,perl=F,fixed=T)
labs <- gsub("(","",labs,perl=F,fixed=T)
labs <- gsub(")","",labs,perl=F,fixed=T)
labs <- strsplit(labs,",")
labs <- round(as.numeric(unlist(labs)),rnd)
nn <- length(labs)/2
labs. <- list()
for(i in 1:nn){
	labs.[[i]] <- labs[c(2*i-1,2*i)]
	labs.[[i]] <- paste("[",paste(labs.[[i]],collapse=", "),")",sep="")
}
labs.[[length(labs.)]] <- gsub(")","]",labs.[[length(labs.)]],fixed=T)
labs. <- unlist(labs.)
return(labs.)}



########################################################
## Prediction diagnostics
########################################################

sepPlot <- function(preds,values){
  predmat <- data.frame(Y=values,FITTED=preds)
  predmat <- predmat[order(predmat$FITTED),]
  predmat <- na.omit(predmat)
  predmat$ORDER <- 1:nrow(predmat)
  par(mar=c(.5,4,.5,.5))
  plot(x=predmat$ORDER,y=predmat$FITTED,col=NA,xaxt="n",ylab="Predicted probabilites")
  psub <- predmat[predmat$Y==1,]
  xs <- jitter(psub$ORDER)
  segments(x0=xs,x1=xs,y0=0,y1=psub$Y,col="grey")
  lines(predmat$ORDER,predmat$FITTED,lwd=2)
}




ROC <- function(preds,values,n=1000){
#   cutz <- seq(min(preds),max(preds),length.out=n)
  cutz <- seq(0,1,length.out=n)
  true.pos <- c()
  false.pos <- c()
  for(i in 1:length(cutz)){
#     tab <- matrix(0,2,2);row.names(tab) <- c(0,1);colnames(tab) <- c(0,1)
#     tab0 <- table(1*(preds<=cutz[i]),values)
#     tab[row.names(tab0)[row.names(tab0)%in%row.names(tab)],] <- tab0[row.names(tab0)[row.names(tab0)%in%row.names(tab)],]
# 
    false.pos[i] <- sum((preds>=cutz[i])*(values==0),na.rm=T)/sum(values==0,na.rm=T)
    true.pos[i] <- sum((preds>=cutz[i])*(values==1),na.rm=T)/sum(values==1,na.rm=T)
#     false.pos[i] <- 1-tab["0","0"]/sum(tab[,"0"])
#     true.pos[i] <- 1-tab["1","1"]/sum(tab[,"1"])
  }
  return(list(false.pos=false.pos,true.pos=true.pos))
}


##########################
## AUC function
##########################

AUC <- function(y, fitted){
  n1 <- sum(y)
  n <- length(y)
  out <- (mean(rank(fitted)[y == 1]) - (n1 + 1)/2)/(n - n1)
  return(out)
}


##########################
## Cutoff
##########################
i<- 1

cutOff <- function(preds,values,n=1000){
  cutz <- seq(min(preds,na.rm=T),max(preds,na.rm=T),length.out=n)
  pntz <- c()
  for(i in 1:length(cutz)){
    tab <- matrix(0,2,2);row.names(tab) <- c("FALSE","TRUE");colnames(tab) <- c(0,1)
    tab0 <- table(preds>cutz[i],values)
    tab[row.names(tab0)[row.names(tab0)%in%row.names(tab)],] <- tab0[row.names(tab0)[row.names(tab0)%in%row.names(tab)],]
    pntz[i] <- sum(diag(tab))/sum(tab)
  }
  return(cutz[pntz==max(pntz,na.rm=T)][1])
}



#################################################
## Nils' functions
#################################################




streg.sort <- function(data, unitvar, timevar) {
  data[order(data[,unitvar], data[,timevar]),]
}

streg.tlag <- function(data, timevar, lagvar, order=1) {
  result <- numeric(nrow(data))
  timevalues <- data[,timevar]
  varvalues <- data[,lagvar]
  times <- unique(timevalues)
  
  # set initial lag values to NA
  for (time in 1:order) {
    result[timevalues==times[time]] <- NA
  }
  
  # compute lag
  for (time in (order+1):length(times)) {
    result[timevalues==times[time]] <- varvalues[timevalues==times[time-order]]
  }
  result
}

streg.slag <- function(data, timevar, lagvar, Wmat) {
  result <- numeric(nrow(data))
  timevalues <- data[,timevar]
  varvalues <- data[,lagvar]
  times <- unique(timevalues)
  for (time in 1:length(times)) {
    result[timevalues==times[time]] <- Wmat %*% varvalues[timevalues==times[time]]
  }
  result
}




##################################
# Prediction functions
##################################


predmanQF <- function(mod,var,vals,quadratic=F,fax=NULL){
  nsim<-10000
  # vars <- unlist(strsplit(as.character(mod$formula)[3],split=" + ",fixed=T))
  vars <- names(coefficients(mod))
  n <- length(vals)
  X <- as.data.frame(matrix(NA,nrow=n,ncol=length(vars)))
  names(X) <- names(coefficients(mod))
  vars. <- vars[!(grepl("as.factor",vars,fixed=T)|grepl("Intercept",vars,fixed=T))]
  X[,vars.] <- as.data.frame(matrix(apply(mod$model[,vars.],2,function(x){median(x,na.rm=T)}),nrow=n,ncol=length(vars.),byrow=T))
  X[,"(Intercept)"] <- 1
  X[,var] <- vals
  if(quadratic==T){
    quadz <- names(mod$model)[grep("I(",names(mod$model),fixed=T)]
    for(k in 1:length(quadz)){
      quadd <- names(mod$model)[apply(as.data.frame(names(mod$model)),1,function(x){grepl(x,quadz[k])})]
      X[,quadz[k]] <- X[,quadd]^2}}
  if(length(fax)>0){
    for(k in 1:length(fax)){
      X[,vars==paste("as.factor(",names(fax[k]),")",fax[[k]],sep="")]<-0
      X[,which(grepl("as.factor",vars)&grepl(names(fax[k]),vars)&(!grepl(fax[[k]],vars)))] <- 0
      X[,which(grepl("as.factor",vars)&grepl(names(fax[k]),vars)&grepl(fax[[k]],vars))] <- 1
    }}
  betas <- rmvnorm(n=nsim,mean=coefficients(mod),sigma=vcov(mod))
  mu <- betas%*%t(X)
  pred.y <- mod$family$linkinv(mu)
  y <- apply(pred.y,2,median)
  y.lo <- apply(pred.y,2,function(x){quantile(x,.025)})
  y.hi <- apply(pred.y,2,function(x){quantile(x,.975)})
  rrs <- 100*(pred.y[,ncol(pred.y)]/pred.y[,1]-1)
  rr <- mean(rrs)
  rr.lo <- quantile(rrs,.025)
  rr.up <- quantile(rrs,.975)
  return(list(y=y,lo=y.lo,up=y.hi, rr=rr,rr.lo=rr.lo,rr.up=rr.up))}

#########################################################

# var <- "CONVOYS_IV"
# xlab <- "X"
# ylab <- "Y"
# vals <- xs
# mod
# var="CONVOYS_VIL_IV"
# xlab="Convoys derailed by partisans"
# ylab="Houses destroyed by Germans"
# cax=.9
# head(mod$model)
# vars.%in%names(mod$model)

require(mvtnorm)
predmanME <- function(mod,var,vals,quadratic=F,fax=NULL){
  nsim<-10000
  # vars <- unlist(strsplit(as.character(mod$formula)[3],split=" + ",fixed=T))
  # vars <- names(mod$model)
  vars <- names(coefficients(mod))
  vars <- gsub(")vec",").vec",vars)
  n <- length(vals)
  X <- as.data.frame(matrix(NA,nrow=n,ncol=length(vars)))
  names(X) <- names(coefficients(mod))
  # names(X) <- vars
  vars. <- vars[!(grepl("as.factor",vars,fixed=T)|(grepl("fitted",vars,fixed=T))|grepl("Intercept",vars,fixed=T))]
  X[,vars.] <- as.data.frame(matrix(apply(mod$model[,vars.],2,function(x){median(x,na.rm=T)}),nrow=n,ncol=length(vars.),byrow=T))
  X[,"(Intercept)"] <- 1
  X[,var] <- vals
  if(quadratic==T){
    quadz <- names(mod$model)[grep("I(",names(mod$model),fixed=T)]
    for(k in 1:length(quadz)){
      quadd <- names(mod$model)[apply(as.data.frame(names(mod$model)),1,function(x){grepl(x,quadz[k])})]
      X[,quadz[k]] <- X[,quadd]^2}}
  if(length(fax)>0){
    for(k in 1:length(fax)){
      X[,vars==paste("as.factor(",names(fax[k]),")",fax[[k]],sep="")]<-0
      X[,which(grepl("as.factor",vars)&grepl(names(fax[k]),vars)&(!grepl(fax[[k]],vars)))] <- 0
      X[,which(grepl("as.factor",vars)&grepl(names(fax[k]),vars)&grepl(fax[[k]],vars))] <- 1
    }}
  for(k in 1:length(grep("fitted",names(X)))){
    X[,grep("fitted",names(X))[k]]<-apply(mod$model[,grep("fitted",names(mod$model))],2,median)[k]}
  betas <- rmvnorm(n=nsim,mean=coefficients(mod),sigma=vcov(mod))
  mu <- betas%*%t(X)
  pred.y <- mod$family$linkinv(mu)
  y <- apply(pred.y,2,median)
  y.lo <- apply(pred.y,2,function(x){quantile(x,.025)})
  y.hi <- apply(pred.y,2,function(x){quantile(x,.975)})
  rrs <- 100*(pred.y[,ncol(pred.y)]/pred.y[,1]-1)
  rr <- mean(rrs)
  rr.lo <- quantile(rrs,.025)
  rr.up <- quantile(rrs,.975)
  return(list(y=y,lo=y.lo,up=y.hi, rr=rr,rr.lo=rr.lo,rr.up=rr.up))}




#########################################################



predmanGAM <- function(mod,var,vals,quadratic=F,fax=NULL){
  nsim<-10000
  vars <- names(mod$model)
  n <- length(vals)
  X <- as.data.frame(matrix(NA,nrow=n,ncol=length(vars)))
  names(X) <- vars
  vars. <- vars[!(grepl("as.factor",vars,fixed=T)|grepl("Intercept",vars,fixed=T))]
  X[,vars.] <- as.data.frame(matrix(apply(mod$model[,vars.],2,function(x){mean(x,na.rm=T)}),nrow=n,ncol=length(vars.),byrow=T))
  # X[,"(Intercept)"] <- 1
  X[,var] <- vals
  for(j in 1:length(vars.)){X[,vars.[j]] <- as.numeric(as.character(X[,vars.[j]]))}
  if(quadratic==T){
    quadz <- names(mod$model)[grep("I(",names(mod$model),fixed=T)]
    for(k in 1:length(quadz)){
      quadd <- names(mod$model)[apply(as.data.frame(names(mod$model)),1,function(x){grepl(x,quadz[k])})]
      X[,quadz[k]] <- X[,quadd]^2}}
  if(length(fax)>0){
    for(k in 1:length(fax)){
      X[,which(grepl(names(fax[k]),vars))] <- as.factor(fax[[k]])
    }}
  X <- as.data.frame(X)
  preds <- predict.gam(mod,newdata=X,type="link",se.fit=T)
  preds <- cbind(preds$fit,preds$fit-1.96*preds$se.fit,preds$fit+1.96*preds$se.fit)
  preds <- mod$family$linkinv(preds)    
  y <- preds[,1]
  y.lo <- preds[,2]
  y.hi <- preds[,3]
  rr <- 100*(y[length(y)]/y[1]-1)
  rr.lo <- 100*(y.lo[length(y.lo)]/y.lo[1]-1)
  rr.up <- 100*(y.hi[length(y.hi)]/y.hi[1]-1)
  return(list(y=y,lo=y.lo,up=y.hi, rr=rr,rr.lo=rr.lo,rr.up=rr.up))}


#####################################################

# var <- "CONVOYS_VIL_IV"
# xlab <- "X"
# ylab <- "Y"
# mod
# var="CONVOYS_IV"
# xlab="Convoys derailed by partisans"
# ylab="Houses destroyed by Germans"
# cax=.9

contplot <- function(mod,var,xlab="X",ylab="Y",main="",cax=1,fax=NULL){
  xs <- seq(min(mod$model[,var],na.rm=T),max(mod$model[,var],na.rm=T),length.out=100)
  preds <- predmanQF(mod,var,vals=xs,fax=fax)
  pred.mean <- preds$y
  pred.lo <- preds$lo
  pred.hi <- preds$up
  par(mar=c(4,4,.5,.5))
  plot(x=xs,y=pred.mean,type="l",ylim=c(min(pred.lo),max(pred.hi)),main=main,ylab=ylab,xlab=xlab,cex.axis=.9,cex.lab=cax,cex.main=1.2,bty="l")
  lines(x=xs,y=pred.lo,lty=2)
  lines(x=xs,y=pred.hi,lty=2)
  rug(jitter(mod$model[,var]))
}



contplot <- function(mod,var,xlab="X",ylab="Y",main="",cax=1,fax=NULL){
  xs <- seq(min(mod$model[,var],na.rm=T),max(mod$model[,var],na.rm=T),length.out=100)
  preds <- predmanQF(mod,var,vals=xs,fax=fax)
  pred.mean <- preds$y
  pred.lo <- preds$lo
  pred.hi <- preds$up
  par(mar=c(4,4,.5,.5))
  plot(x=xs,y=pred.mean,type="l",ylim=c(min(pred.lo),max(pred.hi)),main=main,ylab=ylab,xlab=xlab,cex.axis=.9,cex.lab=cax,cex.main=1.2,bty="l")
  # lines(x=xs,y=pred.lo,lty=2)
  # lines(x=xs,y=pred.hi,lty=2)
  polygon(x=c(xs,rev(xs)),y=c(pred.lo,rev(pred.hi)),col="gray",border="NA")
  lines(x=xs,y=pred.mean)
  rug(jitter(mod$model[,var]))
}




contplot_GAM <- function(mod,var,xlab="X",ylab="Y",main="",cax=1,fax=NULL,quadratic=F,ylim=NULL,col="gray"){
  xs <- seq(min(mod$model[,var],na.rm=T),max(mod$model[,var],na.rm=T),length.out=100)
  preds <- predmanGAM(mod,var,vals=xs,fax=fax,quadratic=quadratic)
  pred.mean <- preds$y
  pred.lo <- preds$lo
  pred.hi <- preds$up
  if(length(ylim)==0){ylim <- c(min(pred.lo),max(pred.hi))}
  par(mar=c(4,4,.5,.5))
  plot(x=xs,y=pred.mean,type="l",ylim=ylim,main=main,ylab=ylab,xlab=xlab,cex.axis=.9,cex.lab=cax,cex.main=1.2,bty="l")
  # lines(x=xs,y=pred.lo,lty=2)
  # lines(x=xs,y=pred.hi,lty=2)
  polygon(x=c(xs,rev(xs)),y=c(pred.lo,rev(pred.hi)),col=col,border="NA")
  lines(x=xs,y=pred.mean)
  #rug(jitter(mod$model[,var]))
}

contplot_ME <- function(mod,var,xlab="X",ylab="Y",main="",cax=1,fax=NULL,quadratic=F){
  xs <- seq(min(mod$model[,var],na.rm=T),max(mod$model[,var],na.rm=T),length.out=100)
  preds <- predmanME(mod,var,vals=xs,fax=fax,quadratic=quadratic)
  pred.mean <- preds$y
  pred.lo <- preds$lo
  pred.hi <- preds$up
  par(mar=c(4,4,.5,.5))
  plot(x=xs,y=pred.mean,type="l",ylim=c(min(pred.lo),max(pred.hi)),main=main,ylab=ylab,xlab=xlab,cex.axis=.9,cex.lab=cax,cex.main=1.2,bty="l")
  # lines(x=xs,y=pred.lo,lty=2)
  # lines(x=xs,y=pred.hi,lty=2)
  polygon(x=c(xs,rev(xs)),y=c(pred.lo,rev(pred.hi)),col="gray",border="NA")
  lines(x=xs,y=pred.mean)
  rug(jitter(mod$model[,var]))
}




#################################################


binplot <- function(mod,var,labname,main,cax=1){
  X <- apply(prematch[,unlist(strsplit(gsub(")","",gsub("s(","",mod$formula,fixed=T))[3],split=" + ",fixed=T))],2,function(x){mean(x,na.rm=T)})
  X <- as.data.frame(matrix(X,nrow=2,ncol=length(X),byrow=T))
  names(X) <- unlist(strsplit(gsub(")","",gsub("s(","",mod$formula,fixed=T))[3],split=" + ",fixed=T))
  X[,var] <- 0:1
  if("DIST_TO_RAIL"%in%names(X)){X$RAIL2 <- X$DIST_TO_RAIL*X$DIST_TO_RAIL}
  
  preds <- predict(mod,newdata=X,type="link",se.fit=T)
  pred.mean <- preds$fit
  pred.lo <- preds$fit-1.96*preds$se.fit
  pred.hi <- preds$fit+1.96*preds$se.fit
  pred.mean <- 1/(1+exp(-pred.mean))
  pred.lo <- 1/(1+exp(-pred.lo))
  pred.hi <- 1/(1+exp(-pred.hi))
  
  par(mar=c(4,4,2,.5))
  plot(x=X[,var],y=pred.mean,type="p",ylim=c(min(pred.lo),max(pred.hi)),main=main,ylab="Prob. of Resettlement",xlab=labname,cex.axis=.9,cex.lab=cax,cex.main=1.2,xlim=c(-2.5,3.5),xaxt="n",pch=1,col=NA)
  axis(1,at=0:1)
  
  segments(x0=seq(X[1,var]-.3,X[1,var]+.3,length.out=100),x1=seq(X[1,var]-.3,X[1,var]+.3,length.out=100),y0=pred.lo[1],y1=pred.hi[1],lty=1,col="grey")
  segments(x0=seq(X[2,var]-.3,X[2,var]+.3,length.out=100),x1=seq(X[2,var]-.3,X[2,var]+.3,length.out=100),y0=pred.lo[2],y1=pred.hi[2],lty=1,col="grey")
  
  segments(x0=X[1,var]-.3,x1=X[1,var]+.3,y0=pred.mean[1],y1=pred.mean[1],lty=1)
  segments(x0=X[2,var]-.3,x1=X[2,var]+.3,y0=pred.mean[2],y1=pred.mean[2],lty=1)
  segments(x0=X[1,var]-.3,x1=X[1,var]+.3,y0=pred.lo[1],y1=pred.lo[1],lty=2)
  segments(x0=X[2,var]-.3,x1=X[2,var]+.3,y0=pred.lo[2],y1=pred.lo[2],lty=2)
  segments(x0=X[1,var]-.3,x1=X[1,var]+.3,y0=pred.hi[1],y1=pred.hi[1],lty=2)
  segments(x0=X[2,var]-.3,x1=X[2,var]+.3,y0=pred.hi[2],y1=pred.hi[2],lty=2)
  CurlyBraces(1.5, (pred.mean[2]-pred.mean[1])/2+pred.mean[1], pred.mean[2]-pred.mean[1], pos = 1, direction = 1 )
  #text(2,(pred.mean[2]-pred.mean[1])/2+pred.mean[1],paste(ifelse(pred.mean[2]-pred.mean[1]>0,"+ ","- "),round(pred.mean[2]-pred.mean[1],2),"\n","(",round(min(c(pred.lo[2]-pred.lo[1],pred.hi[2]-pred.hi[1])),2),",",round(max(c(pred.lo[2]-pred.lo[1],pred.hi[2]-pred.hi[1])),2),")",sep=""),pos=4,cex=.8)
  text(2,(pred.mean[2]-pred.mean[1])/2+pred.mean[1],paste(if(pred.mean[2]-pred.mean[1]>0){"+ "},round(pred.mean[2]-pred.mean[1],2),sep=""),pos=4,cex=.8)
}





###############################
## Correlation matrix legend
###############################

legend.v2 <- function (x, y = NULL, legend, fill = NULL, col = par("col"), border = "black", lty, lwd, pch, angle = 45, density = NULL, bty = "o", bg = par("bg"), box.lwd = par("lwd"), box.lty = par("lty"), box.col = par("fg"), pt.bg = NA, cex = 1, pt.cex = cex, pt.lwd = lwd, xjust = 0, yjust = 1, x.intersp = 1, y.intersp = 1, adj = c(0, 0.5), text.width = NULL, text.col = par("col"), merge = do.lines && has.pch, trace = FALSE, plot = TRUE, ncol = 1, horiz = FALSE, title = NULL, inset = 0, xpd, title.col = text.col, title.adj = 0.5,seg.len = 2){
if (missing(legend) && !missing(y) && (is.character(y) || 
is.expression(y))) {
legend <- y
y <- NULL
}
mfill <- !missing(fill) || !missing(density)
if (!missing(xpd)) {
op <- par("xpd")
on.exit(par(xpd = op))
par(xpd = xpd)
}
title <- as.graphicsAnnot(title)
if (length(title) > 1) 
stop("invalid title")
legend <- as.graphicsAnnot(legend)
n.leg <- if (is.call(legend)) 
1
else length(legend)
if (n.leg == 0) 
stop("'legend' is of length 0")
auto <- if (is.character(x)) 
match.arg(x, c("bottomright", "bottom", "bottomleft", 
"left", "topleft", "top", "topright", "right", "center"))
else NA
if (is.na(auto)) {
xy <- xy.coords(x, y)
x <- xy$x
y <- xy$y
nx <- length(x)
if (nx < 1 || nx > 2) 
stop("invalid coordinate lengths")
}
else nx <- 0
xlog <- par("xlog")
ylog <- par("ylog")
rect2 <- function(left, top, dx, dy, density = NULL, angle, 
...) {
r <- left + dx
if (xlog) {
left <- 10^left
r <- 10^r
}
b <- top - dy
if (ylog) {
top <- 10^top
b <- 10^b
}
rect(left, top, r, b, angle = angle, density = density, 
...)
}
segments2 <- function(x1, y1, dx, dy, ...) {
x2 <- x1 + dx
if (xlog) {
x1 <- 10^x1
x2 <- 10^x2
}
y2 <- y1 + dy
if (ylog) {
y1 <- 10^y1
y2 <- 10^y2
}
segments(x1, y1, x2, y2, ...)
}
points2 <- function(x, y, ...) {
if (xlog) 
x <- 10^x
if (ylog) 
y <- 10^y
points(x, y, ...)
}
text2 <- function(x, y, ...) {
if (xlog) 
x <- 10^x
if (ylog) 
y <- 10^y
text(x, y, ...)
}
if (trace) 
catn <- function(...) do.call("cat", c(lapply(list(...), 
formatC), list("\n")))
cin <- par("cin")
Cex <- cex * par("cex")
if (is.null(text.width)) 
text.width <- max(abs(strwidth(legend, units = "user", 
cex = cex)))
else if (!is.numeric(text.width) || text.width < 0) 
stop("'text.width' must be numeric, >= 0")
xc <- Cex * xinch(cin[1L], warn.log = FALSE)
yc <- Cex * yinch(cin[2L], warn.log = FALSE)
if (xc < 0) 
text.width <- -text.width
xchar <- xc
xextra <- 0
yextra <- yc * (y.intersp - 1)
ymax <- yc * max(1, strheight(legend, units = "user", cex = cex)/yc)
ychar <- yextra + ymax
if (trace) 
catn("  xchar=", xchar, "; (yextra,ychar)=", c(yextra, 
ychar))
if (mfill) {
xbox <- xc * 0.8
ybox <- yc * 0.5
dx.fill <- xbox
}
do.lines <- (!missing(lty) && (is.character(lty) || any(lty > 
0))) || !missing(lwd)
n.legpercol <- if (horiz) {
if (ncol != 1) 
warning("horizontal specification overrides: Number of columns := ", 
n.leg)
ncol <- n.leg
1
}
else ceiling(n.leg/ncol)
has.pch <- !missing(pch) && length(pch) > 0
if (do.lines) {
x.off <- if (merge) 
-0.7
else 0
}
else if (merge) 
warning("'merge = TRUE' has no effect when no line segments are drawn")
if (has.pch) {
if (is.character(pch) && !is.na(pch[1L]) && nchar(pch[1L], 
type = "c") > 1) {
if (length(pch) > 1) 
warning("not using pch[2..] since pch[1L] has multiple chars")
np <- nchar(pch[1L], type = "c")
pch <- substr(rep.int(pch[1L], np), 1L:np, 1L:np)
}
}
if (is.na(auto)) {
if (xlog) 
x <- log10(x)
if (ylog) 
y <- log10(y)
}
if (nx == 2) {
x <- sort(x)
y <- sort(y)
left <- x[1L]
top <- y[2L]
w <- diff(x)
h <- diff(y)
w0 <- w/ncol
x <- mean(x)
y <- mean(y)
if (missing(xjust)) 
xjust <- 0.5
if (missing(yjust)) 
yjust <- 0.5
}
else {
h <- (n.legpercol + (!is.null(title))) * ychar + yc
w0 <- text.width + (x.intersp + 1) * xchar
if (mfill) 
w0 <- w0 + dx.fill
if (do.lines) 
w0 <- w0 + (seg.len + +x.off) * xchar
w <- ncol * w0 + 0.5 * xchar
if (!is.null(title) && (abs(tw <- strwidth(title, units = "user", 
cex = cex) + 0.5 * xchar)) > abs(w)) {
xextra <- (tw - w)/2
w <- tw
}
if (is.na(auto)) {
left <- x - xjust * w
top <- y + (1 - yjust) * h
}
else {
usr <- par("usr")
inset <- rep(inset, length.out = 2)
insetx <- inset[1L] * (usr[2L] - usr[1L])
left <- switch(auto, bottomright = , topright = , 
right = usr[2L] - w - insetx, bottomleft = , 
left = , topleft = usr[1L] + insetx, bottom = , 
top = , center = (usr[1L] + usr[2L] - w)/2)
insety <- inset[2L] * (usr[4L] - usr[3L])
top <- switch(auto, bottomright = , bottom = , bottomleft = usr[3L] + 
h + insety, topleft = , top = , topright = usr[4L] - 
insety, left = , right = , center = (usr[3L] + 
usr[4L] + h)/2)
}
}
if (plot && bty != "n") {
if (trace) 
catn("  rect2(", left, ",", top, ", w=", w, ", h=", 
h, ", ...)", sep = "")
rect2(left, top, dx = w, dy = h, col = bg, density = NULL, 
lwd = box.lwd, lty = box.lty, border = box.col)
}
xt <- left + xchar + xextra + (w0 * rep.int(0:(ncol - 1), 
rep.int(n.legpercol, ncol)))[1L:n.leg]
yt <- top - 0.5 * yextra - ymax - (rep.int(1L:n.legpercol, 
ncol)[1L:n.leg] - 1 + (!is.null(title))) * ychar
if (mfill) {
if (plot) {
fill <- rep(fill, length.out = n.leg)
rect2(left = xt+.22, top = yt + ybox/2, dx = xbox, dy = ybox, 
col = fill, density = density, angle = angle, 
border = border)
}
xt <- xt + dx.fill
}
if (plot && (has.pch || do.lines)) 
col <- rep(col, length.out = n.leg)
if (missing(lwd)) 
lwd <- par("lwd")
if (do.lines) {
if (missing(lty)) 
lty <- 1
lty <- rep(lty, length.out = n.leg)
lwd <- rep(lwd, length.out = n.leg)
ok.l <- !is.na(lty) & (is.character(lty) | lty > 0)
if (trace) 
catn("  segments2(", xt[ok.l] + x.off * xchar, ",", 
yt[ok.l], ", dx=", seg.len * xchar, ", dy=0, ...)")
if (plot) 
segments2(xt[ok.l] + x.off * xchar, yt[ok.l], dx = seg.len * 
xchar, dy = 0, lty = lty[ok.l], lwd = lwd[ok.l], 
col = col[ok.l])
xt <- xt + (seg.len + x.off) * xchar
}
if (has.pch) {
pch <- rep(pch, length.out = n.leg)
pt.bg <- rep(pt.bg, length.out = n.leg)
pt.cex <- rep(pt.cex, length.out = n.leg)
pt.lwd <- rep(pt.lwd, length.out = n.leg)
ok <- !is.na(pch) & (is.character(pch) | pch >= 0)
x1 <- (if (merge && do.lines) 
xt - (seg.len/2) * xchar
else xt)[ok]
y1 <- yt[ok]
if (trace) 
catn("  points2(", x1, ",", y1, ", pch=", pch[ok], 
", ...)")
if (plot) 
points2(x1, y1, pch = pch[ok], col = col[ok], cex = pt.cex[ok], 
bg = pt.bg[ok], lwd = pt.lwd[ok])
}
xt <- xt + x.intersp * xchar
if (plot) {
if (!is.null(title)) 
text2(left + w * title.adj, top - ymax, labels = title, 
adj = c(title.adj, 0), cex = cex, col = title.col)
text2(xt, yt, labels = legend, adj = adj, cex = cex, 
col = text.col)
}
invisible(list(rect = list(w = w, h = h, left = left, top = top), 
text = list(x = xt, y = yt)))
}


###############################
## Correlation matrix plot
###############################


cormat <- function(data,vars,labels=vars){

Xs.main <- data[,vars]
Xs.main <- na.omit(Xs.main)

mat.cor <- round(cor(Xs.main),2)
rownames(mat.cor) <- labels
colnames(mat.cor) <- rownames(mat.cor) 

# Color palette
pal <- c("darkblue","blue","lightblue","lightblue1","grey95","beige","lightgoldenrod1","orange","red","darkred")
cols.mat <- matrix(NA,nrow=nrow(mat.cor),ncol=ncol(mat.cor))
cols.mat <- ifelse(mat.cor>=-1&mat.cor<(-.8),pal[1],
ifelse(mat.cor>=-(.8)&mat.cor<(-.6),pal[2],
ifelse(mat.cor>=-(.6)&mat.cor<(-.4),pal[3],
ifelse(mat.cor>=-(.4)&mat.cor<(-.2),pal[4],
ifelse(mat.cor>=-(.2)&mat.cor<(0),pal[5],
ifelse(mat.cor>=-(0)&mat.cor<(.2),pal[6],
ifelse(mat.cor>=-(.2)&mat.cor<(.4),pal[7],
ifelse(mat.cor>=-(.4)&mat.cor<(.6),pal[8],
ifelse(mat.cor>=-(.6)&mat.cor<(.8),pal[9],
pal[10]
)))))))))
cols.mat <- cols.mat[ nrow(cols.mat):1,]

# pch matrix
pchs.mat <- matrix(NA,nrow=nrow(mat.cor),ncol=ncol(mat.cor))
pchs.mat <- ifelse(mat.cor<0,4,NA)
pchs.mat <- pchs.mat[ nrow(pchs.mat):1,]


par(mar=c(0,0,0,0))
plot(c(-8,nrow(mat.cor)+1), c(0,nrow(mat.cor)+3), type = "n", xlab="", ylab="", asp = 1,axes=F)

ups <- upper.tri(cols.mat,diag=T)[ncol(cols.mat):1,]
cols.mat[which(ups)] <- NA
ups <- upper.tri(pchs.mat,diag=T)[ncol(pchs.mat):1,]
pchs.mat[which(ups)] <- NA

for(i in 1:nrow(cols.mat)){
for(j in 1:ncol(cols.mat)){
points(x=i,y=j,cex=2.4,col=cols.mat[j,i],pch=15)
points(x=i,y=j,cex=2.4,col="grey",pch=pchs.mat[j,i])
}
}
text(x=0.5,y=ncol(cols.mat):1,pos=2,labels = c(NA,labels[-1]),cex=.7)
text(x=1:(ncol(cols.mat))-.35,y=ncol(cols.mat):1,pos=4,labels = c(labels[-length(labels)],NA),cex=.75, srt = 90)

legend.v2(x="bottom",cex=.8,fill=pal,legend=c("[-1, -.8)","[-.8, -.6)","[-.6, -.4)","[-.4, -.2)","[-.2, 0)","[0, .2)","[.2, .4)","[.4, .6)","[.6, .8)","[.8, 1]"),title="Correlation",ncol=5,pch=c(4,4,4,4,4,NA,NA,NA,NA,NA),bty="n",col=c(rep("darkgrey",5),rep(NA,5)))



}



################################################################
## Summary stats
################################################################

# data <- data
# varz <- c("Pub","HARDNEWS_WS")


sumstat <- function(data,varz,labs=NULL){
  sum.mat <- data.frame(Obs=apply(data[,varz],2,function(x){sum(!is.na(x))}),
                        Mean=apply(data[,varz],2,function(x){mean(x,na.rm=T)}),
                        Median=apply(data[,varz],2,function(x){median(x,na.rm=T)}),
                        SD=apply(data[,varz],2,function(x){sd(x,na.rm=T)}),
                        Min=apply(data[,varz],2,function(x){min(x,na.rm=T)}),
                        Max=apply(data[,varz],2,function(x){max(x,na.rm=T)}))
  if(length(labs)==0){labs <- varz}
  row.names(sum.mat) <- labs
  return(sum.mat)
}